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^3 ' Abstract 
^ ■ 

The calculation of the density matrix for fermions and bosons in the Grand 
Q>^ ■ Canonical Ensemble allows an efficient way for the inclusion of fermionic and 

■ bosonic statistics at all temperatures. It is shown that in a Path Integral For- 

o 

OO 



mulation fermionic density matrix can be expressed via an integration over 
a novel representation of the universal temperature dependent functional. 
While several representations for the universal functional have already been 
developed, they are usually presented in a form inconvenient for computer 



o 

CJ ■ calculations. In this work we discuss a new representation for the universal 



• ^ . functional in terms of Hankel functions which is advantageous for computa- 

X 

. tional applications. Temperature scaling for the universal functional and its 

derivatives are also introduced thus allowing an efficient rescaling rather then 
recalculation of the functional at different temperatures. A simple illustration 
of the method of calculation of density profiles in Grand Canonical ensemble 
is presented using a system of noninteracting electrons in a finite confining 
potential. 
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I. INTRODUCTION 



The standard treatment of fermionic and bosonic systems in canonical ensemble using 
Path Integral Formalism is based on the isomorphism between a system of interacting quan- 
tum particles and a set of interacting classical polymer rings [0]. For a system of bosons 
all the ring contributions to the density matrix come with the same sign thus allowing an 
effective Monte-Carlo sampling of the relevant terms [0]. However, for the case of fermions 
situation is more complicated as different ring contributions to the density matrix come with 
alternating signs thus making statistical fluctuations so big that an unreasonably long av- 
erages are required to provide for the meaningful results. While several promising attempts 
were made to circumvent the sign problem there is still much to be done before this 
problem will be solved. 

However, instead of operating in a canonical ensemble for fermionic and bosonic systems 
it might be more advantageous to work in the Grand Canonical Ensemble. To see why, as 
soon as one knows the normal excitations Ei of the system the grand canonical partition 
function is then 

oo 

= ± ^ log{l ± Zexp{-pEi)) (1) 

i=l 

where Z = exp{[3^). The density matrix can then be expressed as 

p(r,r-;A^) =< ^ ^ ^ K > (2) 

where H is the Hamiltonian. By reformulating the problem in this fashion, one need not 
worry about the explicit antisymmetrization or symmetrization of |r > and \r' > states. The 
problem of finding and p{r,r'] 13, fi) is re-expressed in terms of the normal excitations 
of a system. Thus if in the canonical formulation we had interacting fermionic or bosonic 
particles with a density matrix exp{—l3Ho) and a properly symmetrized wave function, in 
the grand canonical formulation we can in principle reduce the problem to a system of 
non interacting particles described by the density matrix z-'^exp{i3H)±i - '^^^^ appears much 
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simpler than the problem in a canonical formulation but there is one problem. The question 
is how to find a spectrum of normal excitations. For fermionic systems the problem is solved 
in principle in the Kohn Sham (KS) density functional formalism [Q. According to (KS) the 
problem of interacting fermions with a Hamiltonian Hq can be mapped into a problem 
of noninteracting particles described by a modified Hamiltonian H which is, in turn, the 
functional of the density (see the review 0) 

H = -— v' +VUr) + Vnir) + Vxcir) (3) 

where Vionir) is an electron- ion interaction potential, 

V„(r) = f !p:^dV (4) 



is the Hartree potential of the electrons and finally, 

Vxc{r) = -T. 5) 

is a universal exchange correlation electron density functional. 

Thus, a self-consistent algorithm for finding the density matrix in the Grand Canonical 
Formulation would be implemented as following 

a) Use p{r, r; (3) to find a new effective Hamiltonian H 

b) Given H calculate a new p(r, r; (3) from 

p{ry-.P.p) =< r\ ^_,^J^^^^^ \r'> (6) 



c) Update the chemical potential p (which comes via Z = exp{l3ii)) 

d) If self-consistency is not achieved go to a) 

Up to this point the algorithm is rather general and can be implemented in many ways. 
One of the implementations of this algorithm for fermions is originally due to S. Goedecker 
and lately extended by Baer The most important step is to implement part b). So 
that 
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—— —-— -|r' >= V V < r\k >< k\—— — — -\k' >< k'lr' > 



(7) 



where < r|fc > is any convenient basis set. By approximating z-^exp{pH)+i ~ Pol{H;[3) 
using Chebyshev Polynomials Goedecker et. al obtained an efficient way of calculating the 
density matrix using Pol{H] (3) as a propagator. 

An alternative and potentially more beneficial method for implementation of step b) 
is to use the Path Integral approach PJTO[]. In the Path Integral Formalism to evaluate 
the density matrix one would insert P auxiliary states thus reducing the problem to the 
evaluation of < '^|(^^;]j(^(^z]:))±i)^k' >• At high temperature 

< >-< r\exp{~l,{H - ,W >~ (8) 

< r\exp{-^Ttin)\r' > exp{--^V{r')) 

so the conventional use of the Trotter approximation for the splitting of the Hamiltonian 
into kinetic and potential energy is applicable. At low temperatures (] —>■ in the case of 
fermion statistics 

<''l' e.p(/?(^/ >-'"^-^> 

where 6{x) is a step function and the propagator becomes invariant with respect to P. So the 
major difference between using the grand canonical approach in the Path Integral Formu- 
lation is the approximation of the intermediate propagator. As a function of P, propagator 
< ^|( e3;p(/3(H-;j))±i ) 1^^ ^ ^^^y ^^^y estimate at high temperature by factoring out the 
kinetic and potential term. At low temperatures the propagator becomes P-independent 
thus substantially hindering its use in calculations and thus new methods need to be imple- 
mented for the evaluation of this propagator. In this work we present a tractable solution 
valid at all temperatures. 
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II. ALL TEMPERATURE SOLUTION 



The problem of evaluation of the fermionic density matrix at zero temperature allows an 



exact solution in the Path Integral Formulation in terms of Bessel functions [rT| . For the 
case of non-zero temperatures the problem has also been solved with the fermionic density 
matrix represented in integral form [T^. 

Here we will show that there exists a more computationally convenient representation 
for the fermionic and bosonic density matrix in terms of the Hankel functions of the first 
kind. 

To calculate the density matrix in the Grand Canonical Ensemble we will exploit the 
property of a meromorphic function been equal to a summation over its poles and residues 
at those poles. 

1 1 1 

exp{P{H - /x)) ± 1 = ^2 ^ „5oo - m - /i) 

where Wn = 7r(2n + 1) for fermions and Wn = '2'n'n for bosons. The density matrix now can 
be written as 

<r| -— ^ \r'>=± ^ ^ ± y <r\ -\r' > (11) 

In this form it is still difficult to make a transition to a Path Integral Formalism. We 
notice, however, that for Wn> ^ 

1 j-+oo 



iwn — P{H — jj) Jo 

and for — m 



r+oo 

i I dtexp{—tWn — it(]{H — /i)) (12) 
Jo 



-"a 



1 /■+00 

i 



-iWn - [3{H - ^) 



r+oo 

/ dtexp{-tWn + it(3{H - //)) (13) 
Jo 



Thus, evaluation of < r| V > is now substantially simphfied and we can change 

to a Path Integral Formulation. 

1 r+oo 

< r\ — -\r' >= —i / dtexp{—tWn) < r\exp{—itp{H — fi))\r' > (14) 

iWn — f3{H — /i) JO 



Element < r\exp{—itl3{H — fi))\r' > is nothing else but a real time propagator. Its evaluation 
in a Path Integral Formalism is trivial and leads to 

< r\exp{-it(3{H - //))|r' >= limt-^r,^t{^^)^exp{^{'^ + ^)) (15) 



where 



/,2^2 1 , T/.„A P-1 

'2m 2 

1=1 



(2 = f ((r - r<'>)' + (r<'> - r<'+'>)' + (r<''-'> - r^) (17) 



P-2 



i=l 



Finally, for iy„ > 



I 1 I / r+oo-ir, ^-imP.3P . .mil ^^-o^Nx / n 

<r — — - r >= —I / dtexp(—tWn)( — r ^ expu — ^ H —)) 18 

and for —Wn 

I 1 I / ■ /■+°°+''' , , imP ,3P , ..mil fiklt 
< f\ TTTTT T >= ^ / atexp{—twn)( — :-—::) ^ exol— ■?(— H )) (19) 

Using the integral representation of the Hankel's functions |13| 

Z 1 /■+00 (j^ \ 2p' 

and summing over all frequencies, it is easy to show that the expression for the density 
matrix can be then rewritten as following, 

p(r, r'; (3, mu) = + j dv^"" {-)- Re{{^ H^^.^z^)}) (21) 



for fermions 



p(r.r';/3.,) = -^(!:^-/dr-(f)¥fle({Z^f(ij)¥-//« (22) 

2 J l-K p Lp 2 ^ 

-TTi Zo x4P_i„(l) 



+ lim {--(-f)^-^//^^_^(.o)}) 



for bosons, where 
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Zn = {{Kl + i^)Llf2 ■ = 7r(2n + l); = 27rn (23) 

and 



Ll = P((r - r(i))2 + ^V^'^ - r^'^'^)' + (r^^-^^ - r')') (25) 

i=l 

In the case of bosons, a term in the summation corresponding to zero frequency Wq should 
be calculated in a limit wq +0. 

III. INTEGRAL REPRESENTATION 

It is also of interest to show how an integral representation of the universal functional 
appears in this formalism. We will outline the major steps on the example of fermions. We 
notice, that for > equations ( [T^JTB| ) can be written as 



1 /-o 

-I 



iWn — (3{H — fl) J-oo 

and for —Wr, 



1 /•+00 

i 



f dtexp{twn + itp{H - /i)) (26) 



r+co 

/ dtexp{-twn + itPiH - n)) (27) 
Jo 



-iWn - P{H - jj,) 
Correspondent propagators are For Wn > 

< r\-. -y >= -I / dtexp{tw^){-—-)-exp{-i{-^ + 28 



and for —Wr, 



1 , / r+oo+ir, ^jjip 3p mil ^^c^NN / X 

< r\ — -\r >= I / dtexp[—tWn)[ — ^) ^ exp(— «(— ^ H —)) (29) 

—iWn — p[H — jjL) Jo+i-n 27rntp 2ht 2m 

Now summation over frequencies can be implemented. For the case of fermions 
J2 < ^\ -^ _ p(^H - fi) ^^' J ^^^^ exp{tLOn)) < r\exp{it(3{H - /i))|r' > (30) 



n=0 """" I^J ""^ n=0 



1 \ p+oo -1 

V < r|— — — -|r' >= i j dt{ Y] exp{-tu)n)) < r\exp{it/3{H - ij,))\r' > 



n=— oo 



(31) 



Making an explicit summation over frequencies and substituting the above equations 
into (1TT|) we get 



where 



1 exp(-7r^) 
'^1 - exp(-27r^) 



= stgn{t)^ ^ (33) 
It is straightforward to check that in the hmit of — > 00 the above expression coincides 



with that of Ref. Ill 



IV. SCALING OF THE UNIVERSAL FUNCTIONAL 

Above we showed that the density matrix for fermions and bosons can be expressed via 
a universal functional 

p(r, r'; (3, /i) = ±^±^ ± J dr'^TiP, P, L|, Kj.) (34) 

From the explicit form of the universal functional derived above one can easily derive a 
scahng law 

F{(3,P,Ll,K],) = (3-'^F{l,P,Ll(3-'',Klp'^) (35) 

This scaling law can be useful if we consider that the integral over the universal functional 
can be rewritten as 

J dr^^F{(3, P, Ll, Kl) = j dLldKl,P{Ll, Kl^)F{P, P, 1% + K^,) (36) 

where 
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ri'Kl^ _ 1 V{r) + V{r'^ ^"^ 
2m ~ 2 



-^ C''>V> + ^V(r'»)) (37) 

i=l 



and P{Lp, Kp^) is the distribution of pairs of variables {Lp,Kp^). So, in principle, it is 
enough to calculate the universal functional once for say (3 — 1.0 and the value of the 
functional at all other temperatures can be evaluated using the scaling law. The integration 
process can be implemented by first calculating the distribution P{L^p, Kp^) (which actually 
depends upon r and r') and then multiplying this distribution on the values of the universal 
functional calculated on a 2D grid shifted by — along the Kp^ axis. 

The averages and thermodynamical derivatives of the generic operator U (r, r') can now 
be determined from 

< U{r, r') J dvdv'piv, r'; /3, //)[/(r, r') (38) 

and derivatives from 

^<^^'^'') >P ^ [ drdr'^P^^l^^ (39) 
d(3 J d(3 

OjJL J OjJL 

where the thermodynamic derivatives of the density matrix can be calculated using the 
scaling relation for the universal functional giving 



ap(r,r-;/3,/i) _ 1 r , , , , aF(l,P,L|.r^i^|>/?^; 

and 



= ±-i^ / dLUKl,P{Ll, X|>J (42) 

Op 2(5 2 J 

^ dF{l,P,Ll(3-l,Kl(5l) laFO^^IrU^ 3P . 1 

^ Ml ~p di?p -^/[i.p.Lpfi ^,Kp(5^n 
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V. A SIMPLE NUMERICAL STUDY 



To investigate the possible usefulness of this method we have performed a study of a 
simple model system consisting of noninteracting fermions in a confining potential of the 
form 

^(0 = (43) 

Here we assume to work in units m = h = 1. This potential allows 5 discrete levels and 
analytical solution for the \1/ functions. M]. The calculations of the density profiles were 



implemented as follows. First, a universal functional was generated for (3 = 1.0, /i = 0.0, 
and P = 20 (Fig. p. One can see a strong oscillatory nature of the potential signifying the 
fermionic nature of the system under consideration. Next step was to generate the distri- 
bution of -P(lvp, Kp^) which is actually the most time consuming part of the calculations. 
So, to find p{r,r) for each value of r a proper P{Lj>, Kp^) distribution was generated using 
Monte-Carlo sampling. The final step of our calculations was in computing the product of 
the universal potential and the PlLjy, Kp^) distribution. By increasing the chemical po- 
tential one moves the PlLj^, Kp^) distribution profile to the right thus coupling to more 
and more oscillations of the universal functional henceforth introducing more particles in a 
system. Thus calculated density profiles for different values of the chemical potential are 
plotted on Fig. |^. For the sake of comparison density profiles obtained from the exact 
solution of Schroedinger equation are also plotted on Fig. ^. Comparing two figures one 
can notice that even for such a small number of subdivisions P = 20 the overall structure 
of the density profiles is the same in both cases. 



VI. CONCLUSIONS 

Grand Canonical formulation of Path Integral Formalism allows an exact treatment of 
fermionic and bosonic systems at any temperature without encountering a complicated prob- 
lem of antisymmetrization or symmetrization of coordinate eigen states. In this formulation 
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density matrix is expressed as an integral over the universal functional FbosjermiP, P, L^, Kp) 
while particular properties of a system come through the variables Lp, Kp. Computationally 
convenient representation of the universal functional is possible in terms of Hankel functions 
of the first kind so one can, in principle, precalculate the Fhosjermif^, P, Kp) and the 
density itself will be determined by a distribution of the non- universal variables Lp, Kp. 
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FIGURES 

FIG. 1. FfermiP = 1-0, P = 20, L'p,Kpg). One can easily notice strong oscillations due to the 
nature of fermionic statistics. 

FIG. 2. Density profiles for different values of the chemical potential calculated using Grand 
Canonical Path Integral Method. 

FIG. 3. Density profiles for different values of the chemical potential calculated using the 
exact solution of Schroedinger equation. 
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